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A  non-equilibrium  molecular  dynamics  method  is  employed  in  order  to  study  various  phenomena 
accompanying  the  generation  and  propagation  of  shock  waves  in  polyurea  (a  micro-phase  segregated 
elastomer).  Several  recent  studies  reported  in  the  literature  suggested  that  polyurea  has  a  relatively  high 
potential  for  mitigation  of  the  effects  associated  with  blast  and  ballistic  impact.  This  behavior  of  polyurea 
is  believed  to  be  closely  related  to  its  micro-phase  segregated  microstructure  (consisting  of  the  so-called 
“hard  domains”  and  a  soft  matrix)  and  to  different  phenomena/processes  (e.g.  inelastic-deformation  and 
energy-dissipation)  taking  place  at,  or  in  the  vicinity  of,  the  shock  front.  The  findings  obtained  in  the 
present  analysis  are  used  to  help  elucidate  the  molecular-level  character  of  these  phenomena/processes. 
In  addition,  the  analysis  yielded  the  shock  Hugoniot  (i.e.  a  set  of  axial  stress  vs.  density/specific-volume  vs. 
internal  energy  vs.  particle  velocity  vs.  temperature  vs.  shock  speed)  material  states  obtained  in  polyurea 
after  the  passage  of  a  shock  wave.  The  availability  of  a  shock  Hugoniot  is  critical  for  construction  of  a  high 
deformation-rate,  large-strain,  high  pressure  material  models  which  can  be  used  within  a  continuum- 
level  computational  analysis  to  capture  the  response  of  a  polyurea-based  macroscopic  structure  (e.g. 
blast-protection  helmet  suspension  pads)  to  blast/ballistic  impact  loading. 
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1.  Introduction 

Due  to  the  use  of  effective  body  armor,  timely  and  effective  evac¬ 
uation,  and  efficient  and  advanced  wound  care,  military  personnel 
can  today  survive  (previously  fatal)  injuries  caused  by  blasts  and 
bullets/fragments.  However,  even  in  the  absence  of  their  visible 
external  injuries,  the  military  personnel  are  often  found  to  suf¬ 
fer  from  traumatic  brain  injury  (TBI)  after  exposure  to  blast/shock 
waves/loads.  Due  to  its  resulting  high  direct  and  indirect  eco¬ 
nomic  costs  to  society  at  large  (through  lost  earning  potential 
of  the  affected  and  the  burden  of  care  imposed  on  their  fam¬ 
ilies),  TBI  has  become  an  important  societal  problem.  This  has 
led  to  an  urgent  need  for  the  development  of  novel  head  protec¬ 
tion  strategies  requiring  in-depth  understanding  of  relationships 
between  the  morphology/structure  (at  different  length  scales)  of 
the  constituent  materials  and  the  activation/effectiveness  of  dif¬ 
ferent  blast-energy  dissipation/dispersion  mechanisms  [1].  Once 
such  knowledge  is  gained,  new  material  systems  with  superior 
blast-energy  dissipation/dispersion  capabilities  can  be  designed, 
synthesized,  and  their  high  shock-mitigation  potential  validated 
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through  the  use  of  test  coupons  [2].  Lastly,  the  performance  of 
newly  synthesized  and  tested  elastomer-based  shock-mitigation 
efficient  materials,  integrated  (as  strategically  placed  functional 
layers)  into  a  helmet/head  assembly  and  subjected  to  blast 
loads,  can  be  investigated  computationally  and/or  experimentally 

[3.4] . 

Extensive  research  carried  out  by  several  groups  [e.g.  5-7] 
established  polyurea  as  a  material  with  an  unusually  high  poten¬ 
tial  for  blast/ballistic-impact  mitigation.  Consequently,  polyurea 
is  frequently  used  as:  (a)  external  and  internal  wall-sidings  and 
foundation  coating  for  buildings  aimed  at  minimizing  the  degree 
of  structure  fragmentation  and,  in  turn,  minimizing  the  extent 
of  the  associated  collateral  damage  in  the  case  of  a  bomb  blast; 
and  (b)  gun-fire/ballistic  resistant  and  explosion/blast  mitigat¬ 
ing  coating/liner  or  inter-layer  in  blast-resistant  sandwich  panels 
for  military  vehicles  and  structures.  In  our  recent  computational 
work  [3,4],  it  was  demonstrated  that  polyurea  can  also  provide 
blast-mitigation  effects  when  used  as  suspension-pad  material  in 
advanced  combat  helmet  applications.  The  work  presented  in  Refs. 

[3.4]  clearly  identified  a  need  for  a  better  understanding  of  the 
mechanical  response  of  polyurea  under  high  deformation  rate, 
large  strain  high-pressure  loading  conditions  encountered  during 
blast  loading  and  for  a  better  understanding  of  shock-wave  gen¬ 
eration,  propagation  and  structure  as  well  as  for  the  nature  of 
energy-dissipative  processes  at  the  shock  front,  in  this  material. 
In  the  present  work,  an  attempt  is  made  to  get  a  deeper  insight 
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Fig.  1.  A  schematic  of  the  simplified  molecular-level  structure  of  segmented 
polyurea. 


into  the  various  shock  related  phenomena  in  polyurea  via  the  use 
of  molecular-level  computational  analyses. 

Polyurea  is  an  elastomeric  co-polymer  formed  by  the  rapid 
chemical  reaction  between  an  isocyanate  (organic  chemical  con¬ 
taining  isocyanate  -N=C=0  groups)  and  an  amine  (organic 
chemical  containing  amine  -NH2  groups).  Due  to  the  presence  of 
strongly  polar  “urea  linkages”,  polyurea  is  typically  nano-/micro- 
phase  segregated  into  the  so-called  “hard  domains”  (a  high  glass 
transition  temperature,  Tg,  hydrogen-bonded  phase)  and  a  soft 
matrix  (a  low  Tg  compliant  phase).  Since  hydrogen  bonding  pro¬ 
vides  inter-chain  joining,  polyurea  is  often  referred  to  as  being 
a  thermo-plastically  cross-linked  (in  contrast  to  more  commonly 
covalently  cross-linked)  polymer.  A  schematic  is  provided  in  Fig.  1 
in  order  to  identify  the  basic  molecular-level  structural  elements 
in  polyurea. 

As  mentioned  earlier,  the  main  objective  of  the  present  work 
is  to  investigate  various  shock-wave  related  phenomena  using 
molecular-level  computational  methods. 

A  shock  wave  (or  simply  a  shock)  is  a  wave  which  propagates 
through  a  medium  at  a  speed  higher  than  the  sound  speed  and  its 
passage  causes  an  abrupt  and  discontinuous  change  in  the  material 
state  variables  (e.g.,  pressure,  internal  energy,  density,  temperature 
and  particle  velocity).  The  magnitude  of  the  state-variable  changes 
and  the  shock  speed  increase  with  the  strength  of  the  shock.  While 
acoustic  waves  give  rise  to  isentropic  changes  in  the  material  state 
variables,  passage  of  a  shock  is  typically  associated  with  irreversible 
(entropy-increasing)  changes  in  the  same  variables.  The  reason 
behind  this  difference  is  that  shock  involves  very  high  strain  rates 
that  bring  energy-dissipative  viscous  effects  into  prominence. 

A  review  of  the  literature  carried  out  as  a  part  of  the  present 
work  revealed  that  the  molecular-level  computational  methods 
were  first  employed  to  study  shock  waves  more  than  30  years 
ago  [8-1 1  ].  While  the  initial  studies  focused  on  shock  phenomena 
in  dense  fluids,  subsequent  work  also  included  single-crystalline 
solids.  The  key  findings  related  to  the  generation  and  propagation 
of  shocks  in  these  solids  can  be  summarized  as  follows:  (a)  these 


phenomena  are  inherently  more  complex  in  solids  than  in  fluids, 
because  solids,  in  addition  to  the  lattice  parameter,  introduce  a 
new  length  scale  (i.e.  the  size/spacing  of  defects)  which  tends  to 
control  the  nature/extent  of  inelastic-deformation  and  dissipative 
processes  at  the  shock  front;  (b)  shock  propagation  typically  results 
in  the  formation  of  steady  (time-invariant)  waves,  and  this  process 
is  facilitated  by  the  transverse  displacement  of  atoms  which  pro¬ 
duce  inelastic  deformation.  This  deformation  involves  concerted 
slippage  of  atoms  over  each  other  (and  not  viscous  flow  as  in 
the  case  of  shocks  in  fluids);  (c)  in  order  to  eliminate  free-surface 
effects,  molecular-level  modeling  of  shock  is  typically  carried  out 
using  computational  systems  with  periodic  boundary  conditions 
(at  least  in  directions  transverse  to  the  shock- wave  propagation 
direction).  To  attain  steady  wave  conditions  of  the  shock,  compu¬ 
tational  domains  sufficiently  long  in  the  direction  of  shock  wave 
propagation  have  to  be  employed.  Furthermore,  lateral  dimen¬ 
sions  of  the  computational  domain  have  to  be  also  sufficiently 
large  to  avoid  spurious  effects  associated  with  the  use  of  the  peri¬ 
odic  boundary  conditions.  Consequently,  computational  domains 
involving  several  tens  to  hundred  thousands  of  atoms  have  to  be 
employed.  The  corresponding  computational  times  (controlled  by 
the  shock  wave  transversal  time)  are  then  on  the  order  of  5-20  ps. 
This  limits  the  shock  thicknesses  to  around  lOnm  and  the  corre¬ 
sponding  rise  times  to  ca.  0.5-1  ps.  Thus,  weak  shock  waves  with 
thicknesses  of  hundreds  of  nanometers  could  not  be  analyzed  using 
molecular-level  methods  (at  least  in  their  steady-wave  regime). 

While  recognizing  the  aforementioned  aspects  and  potential 
limitations  of  molecular-level  modeling  of  shock,  a  preliminary 
computational  study  of  shock  generation/propagation  in  polyurea 
is  carried  out.  There  are  two  main  objectives  of  the  present  work: 

(a)  To  determine  the  shock  Hugoniot  (centered  on  the  initial  stress- 
free  quiescent  state)  of  polyurea.  A  Hugoniot  is  the  locus  of 
axial  stress  vs.  specific  volume  vs.  energy  density  (vs.  particle 
velocity  vs.  shock  speed)  shocked  (“upstream”)  material  states. 
The  shock  Hugoniot  (or  simply,  the  Hugoniot)  is  often  used  in 
the  derivation  of  the  continuum-level  material  models  (partic¬ 
ularly  in  the  derivation  of  the  equation  of  state)  which  is  used  in 
the  computational  investigations  of  the  response  of  structures 
to  shock  loading.  In  situations  in  which  one  is  interested  only 
in  the  problem  of  planar  shock  propagation/interaction  (in  the 
presence  of  uni-axial  strain  deformation  states),  a  complete  def¬ 
inition  of  the  continuum-level  material  model  is  not  required. 
Instead,  the  knowledge  of  the  corresponding  Hugoniot  relations 
is  sufficient;  and 

(b)  The  Hugoniot  relations  mentioned  in  (a)  provide  a  global 
statement  of  mass,  momentum,  and  energy  conservation 
accompanying  shock-induced  material  transition  from  a  given 
initial  (“downstream”)  equilibrium  state  to  all  possible  final 
(“upstream”)  equilibrium  states  for  steady  planar  shock  waves 
(of  different  strengths).  However,  these  relations  provide  no 
information  about  the  structure  of  the  shock  front  or  the  nature 
of  the  dissipative  structural  rearrangement  mechanisms  that 
lead  to  a  steady  shock  wave.  Hence,  the  second  objective  of 
the  present  work  is  to  carry  out  a  detailed  examination  of 
the  downstream,  shock-front  and  upstream  material  states  (as 
represented  by  the  local  stresses,  strains,  densities/specific  vol¬ 
umes,  temperatures,  etc.)  and  molecular-level  morphology  in 
order  to  identify  and  characterize  these  processes. 

The  organization  of  the  paper  is  as  follows:  a  brief  overview 
of  the  molecular-level  computational  procedure  including  the 
computational  cell  construction,  force  field  identification,  com¬ 
putational  method(s)  selection,  shock-wave  generation  and  the 
problem  definition  are  respectively  presented  in  Sections  2.1 -2.5. 
The  key  results  obtained  in  the  present  work  are  presented  and  dis- 
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cussed  in  Section  3.  A  summary  of  the  work  carried  out  and  the  key 
results/conclusions  are  given  in  Section  4. 

2.  Molecular-level  computational  procedure 

As  mentioned  earlier,  molecular-level  computational  methods 
have  been  employed  in  the  present  work  in  order  to  investigate 
shock  formation  and  propagation  in  polyurea.  Within  these  meth¬ 
ods,  all  atoms  and  bonds  are  explicitly  accounted  for  and  molecular 
mechanics,  dynamics  or  Monte  Carlo  algorithms  are  used  to  quan¬ 
tify  the  behavior  of  the  material  under  investigation. 

While  ab  initio  quantum  mechanics  methods  have  the  advan¬ 
tage  over  the  molecular-level  methods  since  they  do  not  require 
parameterization,  they  have  a  serious  short-coming.  Namely,  due 
to  prohibitively  high  computational  cost,  they  can  be  currently 
employed  only  for  systems  containing  no  more  than  a  few  hundred 
atoms/particles.  As  will  be  shown  below,  while  ab  initio  quantum- 
mechanics  calculations  are  not  directly  used  in  the  present  work, 
some  of  the  computational  ab  initio  quantum  mechanics  results 
are  used  in  the  parameterization  of  the  material  model  at  the 
molecular  length/time  scale.  Utility  of  the  molecular-level  compu¬ 
tational  results  is  greatly  dependent  on  accuracy  and  fidelity  of  the 
employed  force  field  (the  mathematical  expression  which  describe 
various  bonding  and  non-bonding  interaction  forces  between  the 
constituents  of  the  molecular-scale  model).  In  the  present  work, 
the  so-called  “COMPASS”  (Condensed-phase  Optimized  Molecu¬ 
lar  Potentials  for  Atomistic  Simulation  Studies)  force  field  is  used 
[12,13].  This  highly  accurate  force  field  is  of  an  ab  initio  type  since 
most  of  its  parameters  were  determined  by  matching  the  predic¬ 
tions  made  by  the  ab  initio  quantum  mechanics  calculations  to  the 
condensed-matter  experimental  data.  Hence,  it  should  be  recog¬ 
nized  that  the  COMPASS  force  field  is  a  prime  example  of  how  the 
highly  accurate  results  obtained  on  one  length/time  scale  (quan¬ 
tum  mechanic/electronic,  in  the  present  case)  and  the  experimental 
data  can  be  combined  to  parameterize  material  models  used  at 
coarser  length/time  scale  (the  molecular  length/time  scale,  in  the 
present  case). 

Formulation  of  a  molecular-level  simulation  problem  requires, 
at  a  minimum,  specification  of  the  following  three  aspects:  (a)  a 
molecular-level  computational  model  consisting  of  atoms,  ions, 
functional  groups  and/or  molecules;  (b)  a  set  of  force  field  func¬ 
tions;  and  (c)  a  computational  method(s)  to  be  used  in  the 
simulation.  More  details  of  these  three  aspects  of  the  molecular- 
level  modeling  and  simulations  of  polyurea  are  provided  below. 

2.1.  Computational  model 

The  first  step  in  the  molecular-level  computational  analysis  was 
the  construction  of  a  single  polyurea  molecule.  This  was  carried  out 
using  the  Visualizer  [14]  program  from  Accelrys.  A  close-up  view 
of  a  typical  polyurea  molecule  obtained  is  shown  in  Fig.  2(a).  For 
clarity,  different  atomic  species  and  functional  groups  are  labeled 
in  this  figure.  Also,  a  zoomed-out  view  of  a  longer  segment  of  the 
single  polyurea  molecule  is  shown  as  an  inset  in  Fig.  2(a).  It  should 
be  noted  that  the  polyurea  molecule  analyzed  in  the  present  work  is 
based  on  an  isocyanate  consisting  of  di-phenyl  methane  functional 
group  and  a  diamine  consisting  of  poly-tetramethyleneoxide-di- 
phenyl  functional  group  (consisting  of  14  C2H80  units),  to  comply 
with  the  type  of  polyurea  being  investigated  experimentally  in  our 
ongoing  work  [2]. 

In  order  to  model  the  behavior  of  bulk  polyurea  (i.e.,  in  the 
absence  of  any  free-surface  effects)  at  the  molecular  length-scale,  a 
rectangular  box-shaped  computational  unit-cell  is  constructed  and 
the  periodic  boundary  conditions  are  applied  across  the  faces  of  the 
cell.  The  three  edges  (a,  b  and  c)  of  the  cell  are  aligned  respectively 
with  the  three  coordinate  axes  (x,y  and  z).  The  molecular  configura¬ 


tion  within  the  cell  containing  polyurea  chains  whose  construction 
was  described  earlier  is  generated  using  the  following  procedure: 

(a)  The  molecular  chain  constructed  above  is  first  grown  by  a  copy- 
and-attach  process  in  order  to  increase  its  length/molecular 
weight; 

(b)  The  single  polyurea  chain  constructed  in  (a)  is  next  used 
within  the  Amorphous  Cell  program  from  Accelrys  [15]  to  fill 
the  parallelepiped-shaped  computational  cell  of  a  preselected 
size  (7.35  nm  x  1.84  nm  x  1.84  nm)  while  attaining  the  target 
ambient-pressure  density  (1300 kg/m3)  of  polyurea.  The  size 
of  the  computational  cell  was  selected  in  such  a  way  that  the 
following  conditions  are  satisfied:  (i)  In  the  direction  of  shock 
propagation  (direction  x),  the  computational  cell  size  was  set  to 
a  minimum  value  which  is  both  larger  than  the  shock  thickness 
and  sufficiently  large  to  ensure  the  formation  of  a  steady  shock 
wave;  and  (ii)  in  the  lateral  (y  and  z)  directions,  the  size  of  the 
computational  cell  was  set  to  a  minimal  value  which  enables  a 
realistic  analysis  of  the  (material)  transverse  flow;  and 

(c)  The  molecular  configuration  obtained  in  (b)  is  optimized  by 
minimizing  its  potential  energy  with  respect  to  the  positions 
of  the  constituent  atoms. 

An  example  of  the  molecular-level  computational-cell  used  in 
the  present  work  is  displayed  in  Fig.  2(b).  This  cell  contains  1120 
atoms  of  carbon,  1104  atoms  of  hydrogen,  128  atoms  of  nitrogen 
and  192  atoms  of  oxygen. 

As  discussed  earlier,  as-synthesized  polyurea  is  typically  micro¬ 
phase  segregated  and  consists  of  hard  domains  and  a  soft  matrix. 
In  contrast,  the  atomic  scale  polyurea  model  as  depicted  in  Fig.  2(b) 
represents  a  completely  mixed  state  of  this  material.  As  discussed 
in  our  prior  work  [1  ],  within  the  atomic  scale  computational  meth¬ 
ods  used,  micro-phase  separation  processes  cannot  be  generally 
investigated.  The  main  reasons  for  this  are  as  follows: 

(a)  Within  the  molecular  statics,  local  minimization  algorithms  are 
employed  which  are  capable  of  only  finding  a  minimum  energy 
configuration  of  the  given  state  of  the  material  and  typically  do 
not  involve  large-scale  displacement  of  the  chain  segments  or 
chains;  and 

(b)  Molecular  dynamics  is  dominated  by  high  frequency  atomic 
vibrations  and  unrealizable  computational  times  are  needed  to 
simulate  phase  separation  using  this  approach. 

Hence,  atomic-scale  calculations  based  on  the  unit  cell  dis¬ 
played  in  Fig.  2(b)  can  only  provide  information  about  the  fully 
mixed  material  state/properties,  i.e.,  information  about  the  mate¬ 
rial  which  resembles  the  soft-matrix  phase  of  polyurea.  It  should 
be  recalled  that,  within  the  soft  matrix,  hard  and  soft  segments  are 
well  mixed  and,  hence,  the  atomic  scale  material  model  depicted 
in  Fig.  2(b)  is  a  fairly  realistic  representation  of  this  phase  of  the 
polyurea. 

To  gain  some  insight  into  the  molecular-level  response  of 
micro-phase  segregated  polyurea,  a  procedure  is  developed  for 
rearrangement  of  the  polyurea  topology  so  that  the  average  dis¬ 
tance  between  the  urea  linkages  is  greatly  reduced.  This,  in  turn, 
resulted  in  a  substantially  larger  contribution  of  hydrogen  bonding 
to  the  non-bond  potential  energy  of  polyurea.  An  example  of  a  typ¬ 
ical  local  atomic  arrangement  in  “fully-mixed"  and  the  “ segregated " 
polyurea  is  depicted  in  Fig.  2(c)  and  (d),  respectively.  The  urea  link¬ 
ages  can  be  readily  identified  in  these  figures  by  the  presence  of 
the  blue-colored  nitrogen  atoms.  It  should  be  noted  that  although 
the  polyurea  topology  displayed  in  Fig.  2(c)  and  (d)  is  associated 
with  the  segregation  of  hard  and  soft  segments,  this  segregation 
takes  place  at  the  few-nanometer  length  scale.  This  length-scale  is 
substantially  lower  than  its  counterpart  (tens  of  nanometers)  in  the 
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Fig.  2.  (a)  A  segment  of  a  single  polyurea  molecule  with  the  inset  showing  a  zoomed-out  view  of  a  longer  segment  of  the  same  molecule;  (b)  a  molecular-level  computational 
cell  in  which,  for  clarity,  each  molecule  is  assigned  a  unique  color;  (c)  and  (d)  local-level  molecular  structure  in  the  “fully  mixed ”  and  “ micro-segregated ”  polyurea,  respectively. 


case  of  the  real  polyurea.  These  differences  are  expected  to  affect 
shock-generation  and  propagation  processes  in  polyurea  and  will 
be  investigated  in  our  future  work. 

2.2.  Force-fields 


E end-bond-torsion *  stretch-torsion  interactions  between  a  dihedral 
angle  and  its  middle  bond,  Emiddie  bond.torsion,  bend-torsion  inter¬ 
actions  between  a  dihedral  angle  and  one  of  its  valence  angles, 
E angle-torsion*  and  bend-bend-torsion  interactions  between  a  dihe¬ 
dral  angle  and  its  two  valence  angles,  Eangie_angie_torsion ,  terms  as: 


As  stated  above,  the  behavior  of  a  material  system  at  the 
molecular-level  is  governed  by  the  appropriate  force-fields  which 
describe,  in  an  approximate  manner,  the  potential  energy  hyper¬ 
surface  on  which  the  atomic  nuclei  move.  In  other  words,  the 
knowledge  of  force-fields  enables  determination  of  the  potential 
energy  of  a  system  in  a  given  configuration.  In  general,  the  poten¬ 
tial  energy  of  a  system  of  interacting  particles  can  be  expressed  as 
a  sum  of  the  valence  (or  bond),  Evaience,  cross-term,  ECWss-term,  and 
non-bond,  Enon_bond ,  interaction  energies  as: 

Etotal  —  Evaience  +  Ecr0ss-term  +  Enon_bond  (1) 

The  valence  energy  generally  includes  a  bond  stretching  term, 
Ebond*  a  two-bond  angle  term,  Eang[e ,  a  dihedral  bond-torsion  term, 
Etorsion*  an  inversion  (or  an  out-of-plane  interaction)  term,  Eoop, 
and  a  Urey-Bradlay  term  (involves  interactions  between  two  atoms 
bonded  to  a  common  atom),  EUB,  as 


Ecross-term  —  Ebond_bond  +  E  angle-angle  +  Ebond_angie  + 

+^mi'dd/e  _ bond-torsion  +  Eangie_torsion  +  T^ang/e-ang/e-tors/on 


(3) 


The  non-bond  interaction  term,  Enon.bond,  accounts  for  the  inter¬ 
actions  between  non-bonded  atoms  and  includes  the  van  der  Waals 
energy,  Evdw,  and  the  Coulomb  electrostatic  energy,  ECouiomb ,  as: 


Enon-bond  ~  Evdw  +  ^Cou/omb 


(4) 


As  mentioned  earlier,  the  present  molecular-level  analysis  of 
polyurea  employs  the  COMPASS  [12,13]  force-field  for  various  bond 
and  non-bond  interaction  energies  appearing  in  Eqs.  ( 1  )-(4).  A  sum¬ 
mary  of  the  COMPASS  force-field  functions  can  be  found  in  our 
previous  work  [16]. 


2.3.  Computational  method 


mce  ~  Ebond  +  E, 


4-  E+, 


(2) 


The  cross-term  interacting  energy,  ECWss-term,  accounts  for  the 
effects  such  as  bond  length  and  angle  changes  caused  by  the 
surrounding  atoms  and  generally  includes:  stretch-stretch  inter¬ 
actions  between  two  adjacent  bonds,  Ebond_bond,  stretch-bend 
interactions  between  a  two-bond  angle  and  one  of  its  bonds, 
Ebond-angie »  bend-bend  interactions  between  two  valence  angles 
associated  with  a  common  vertex  atom,  Eangie.angie ,  stretch-torsion 
interactions  between  a  dihedral  angle  and  one  of  its  end  bonds, 


Both  molecular  statics  and  molecular  dynamics  simulations 
were  employed  in  the  present  work.  Within  the  molecular  statics 
approach,  the  unit-cell  potential  energy  (as  defined  by  Eqs.  (1  )-(4)) 
is  minimized  with  respect  to  the  position  of  the  constituent  parti¬ 
cles/atoms.  The  potential  energy  minimization  within  Discover  [17] 
(the  atomic  simulation  program  from  Accelrys  used  in  the  present 
work)  is  carried  out  by  combining  the  Steepest  Descent,  Conju¬ 
gate  Gradient  and  the  Newton’s  minimization  algorithms.  These 
algorithms  are  automatically  inactivated/activated  as  the  atomic 
configuration  is  approaching  its  energy  minimum  (i.e.  the  Steepest 
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Descent  method  is  activated  at  the  beginning  of  the  energy  mini¬ 
mization  procedure  while  the  Newtons’  method  is  utilized  in  the 
final  stages  of  this  procedure). 

As  will  be  discussed  in  greater  detail  in  Section  3,  molecular  stat¬ 
ics  is  employed  to  determine  the  state  of  the  material  swept  by  a 
shock.  As  will  be  shown,  this  procedure  is  based  on  the  use  of  (bond¬ 
ing  and  non-bonding)  potential  energy  components  and  neglects 
shock-induced  changes  in  the  (configurational)  entropy  of  the  sys¬ 
tem.  To  assess  the  consequence  of  this  simplification,  the  approach 
described  in  Ref.  [18]  was  considered.  This  approach  defines  a 
dimensionless  parameter  and  states  that  when  this  parameter  is 
significantly  smaller  than  unity,  entropy  effects  can  be  neglected. 
Unfortunately,  detailed  temperature  and  pressure  dependencies  of 
the  material  mechanical  response  of  polyurea  needed  to  evaluate 
this  parameter  were  not  available  with  sufficient  fidelity.  Hence, 
the  results  obtained  by  the  application  of  this  procedure,  which 
suggest  that  the  entropy  effects  may  not  be  highly  critical,  cannot 
be  accepted  with  a  high  level  of  confidence. 

Within  the  molecular  dynamics  approach,  gradient  of  the  poten¬ 
tial  energy  with  respect  to  the  particle  positions  is  first  used  to 
generate  forces  acting  on  the  particles  and,  then,  the  associated 
Newton’s  equations  of  motion  (for  all  particles)  are  integrated 
numerically  in  order  to  track  the  temporal  evolution  of  the  parti¬ 
cle  positions.  Both  the  equilibrium  and  non-equilibrium  molecular 
dynamics  methods  are  employed  in  the  present  work.  Within 
the  equilibrium  molecular-dynamics  methods,  the  system  under 
consideration  is  coupled  to  an  (external)  environment  (e.g.  con¬ 
stant  pressure  piston,  constant  temperature  reservoir,  etc.)  which 
ensures  that  the  system  remains  in  equilibrium  (i.e.  the  system  is 
not  subjected  to  any  thermodynamic  fluxes).  As  will  be  discussed 
in  next  section,  NVT  (where  N  is  the  (fixed)  number  of  atoms,  V , 
the  computational  cell  volume  (also  fixed),  and  T=(298I<)  is  the 
temperature)  equilibrium  molecular  dynamics  is  employed  in  the 
first  stage  of  the  shock  generation  procedure.  In  addition,  as  will 
be  discussed  in  Section  3,  NVE  (. E  is  the  total  energy)  equilibrium 
molecular  dynamics  is  also  employed  during  determination  of  the 
shock  Hugoniot.  Within  non-equilibrium  molecular  dynamics,  the 
system  is  subjected  to  large  perturbations  (finite  changes  in  the 
axial  parameter  of  the  computational  cell,  in  the  present  case) 
which  create  a  thermodynamic  flux  (i.e.  the  flux  of  energy  and 
momentum,  in  the  present  case).  More  details  regarding  the  use 
of  Discover  to  carry  out  molecular  statics  and  molecular  dynamics 
analyses  can  be  found  in  our  prior  work  [1  ]. 

2.4.  Shock-wave  generation 

To  generate  a  planar  shock  (or  more  precisely  a  pair  of  planar 
shocks)  within  the  computational  cell,  the  following  procedure  is 
employed: 

(a)  At  the  beginning  of  the  analysis,  a  “sufficiently  long”  NVT  molec¬ 
ular  dynamics  simulation  is  carried  out  in  order  to  equilibrate 
the  system/material. 

(b)  The  shock  is  then  initiated  (and  driven)  by  continuously  con¬ 
tracting  the  computational  cell  x-direction  lattice  parameter  a 
as: 

a(t)  =  a(t  =  0)  -  2upt  (5) 

where  t  denotes  time,  up  is  the  so-called  “piston”  velocity  (or  equiv¬ 
alently  the  particles  upstream  velocity)  in  the  x-direction.  up  is 
varied  over  a  range  between  187.5  and  1500  m/s  in  order  to  simu¬ 
late  the  generation  and  propagation  of  shock  of  various  strengths. 
Meanwhile,  computational-cell  transverse  lattice  parameters  b  and 
c  are  kept  constant  in  order  to  obtain  planar  (uniaxial-strain)  shock 
conditions.  In  this  process,  the  computational  cell  faces  normal 


Fig.  3.  A  schematic  of  the  generation  of  a  pair  of  shocks  in  a  molecular-level  system 
via  the  process  of  computational-cell  parameter  contraction. 

to  the  shock  propagation-direction  behave  very  similarly  to  the 
impact-surface  of  a  plate-like  target  subjected  to  a  so-called  sym¬ 
metric  “flyer-plate”  impact  test  [19].  The  procedure  employed  here 
generates  a  pair  of  shock  waves  which  propagate,  at  a  shock  speed 
Us,  from  the  cell  boundaries  towards  its  center.  As  schematically 
shown  in  Fig.  3,  these  shock  waves  leave  behind  a  “shocked”  mate¬ 
rial  state  characterized  by  a  higher  material  mass  density  (as  well 
as  by  higher  levels  of  the  internal  energy,  stress,  particle  velocity, 
temperature  and  entropy). 

The  aforementioned  procedure  for  shock-wave  generation  and 
the  subsequent  molecular  statics/dynamics  analyses  are  carried  out 
through  the  use  of  a  Discover  input  file  [17]  which  is  written  in  a 
BTCL  (Basic  Tool  Command  Language)  language.  This  enabled  the 
use  of  a  scripting  engine  that  provides  very  precise  control  of  sim¬ 
ulation  jobs,  e.g.  a  cell  deformation  to  be  carried  out  in  small  steps 
each  followed  by  a  combined  energy-minimization/molecular- 
dynamics  simulation  run. 

2.5.  Problem  formulation 

The  problem  addressed  in  the  present  work  involved  gen¬ 
eration  of  shock  waves  of  different  strengths  (using  the  afore¬ 
mentioned  computational  cell  parameter  contraction  method), 
determination  of  the  associated  shock-Hugoniot  relations  and 
identification  and  elucidation  of  the  main  molecular-level  inelastic- 
deformation/energy  dissipation  processes  taking  place  at  or  in  the 
vicinity  of  the  shock  front.  The  procedure  for  shock  wave  generation 
was  presented  in  the  previous  section. 

As  far  as  the  shock  Hugoniot  determination  is  concerned,  it 
entailed  the  knowledge  of  the  shock-wave  profiles  (and  their 
temporal  evolution)  for  the  axial  stress,  material  density,  particle 
velocity,  internal  energy  and  temperature.  The  latter  are  obtained 
by  lumping  particles/atoms  and  their  (bond  and  non-bond)  poten¬ 
tial  and  kinetic  energy  contributions  into  fixed-width  bins,  in  the 
order  of  their  axial  coordinates.  As  will  be  shown  in  the  next 
section,  two  types  of  bins  are  used:  (a)  a  Lagrangian-type  which 
is  fixed  to  the  initial/reference  state  of  the  computational  cell 
and  (b)  a  moving-type  which  is  attached  to  the  advancing  shock 
front. 

Identification  of  the  molecular-level  inelastic- 
deformation/energy  dissipation  processes  entailed  a  close 
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Fig.  4.  Temporal  evolution  of  the  molecular  level  material  microstructure  accom¬ 
panying  generation  and  propagation  of  a  pair  of  planar  shocks. 


examination  of  the  changes  in  a  material  bond  structure  and 
topology  caused  by  the  passage  of  the  shock. 

3.  Results  and  discussion 

3.1.  Shock-wave  observation  and  structure  characterization 


800 


0  12  3  4 

Distance  in  x-direction,  nm 


An  example  of  the  typical  results,  obtained  in  the  present  work, 
pertaining  to  material  molecular-level  microstructure/topology 
evolution  caused  by  a  continuous  axial  contraction  of  the  com¬ 
putational  cell  is  displayed  in  Fig.  4(a)-(d).  The  results  displayed 
in  these  figures  clearly  reveal  the  generation  of  a  pair  of  planar 
shock  waves  at  the  two  y-z  faces  of  the  computational  cell,  Fig.  4(a), 
and  their  subsequent  propagation  towards  the  center  of  the  cell, 
Fig.  4(b)-(d).  An  approximate  location  of  the  center-point  of  the 
two  shocks  is  indicated  using  arrowheads  in  Fig.  4(a)-(d).  While 
the  two  shocks  collide  at  a  later  simulation  time,  a  shock-wave 
interaction  investigation  is  beyond  the  scope  of  the  present  work. 
The  results  displayed  in  Fig.  4(a)-(d)  show  that  the  shock  waves 
remain  fairly  planar  during  their  motion. 

While  the  results  displayed  in  Fig.  4(a)-(d)  provide  clear  evi¬ 
dence  for  the  formation  and  propagation  of  a  pair  of  opposing 
planar  shock  waves,  they  do  not  offer  any  information  about  the 
structure/shape  of  the  shock-wave  front  or  about  the  state  of  the 
(upstream)  material  swept  by  the  shock.  The  latter  aspects  of  shock- 
wave  generation  and  propagation  within  polyurea  are  addressed  in 
the  remainder  of  this  section  and  in  the  subsequent  sections. 

To  reveal  the  structure/shape  of  the  shock  wave,  the  method  of 
(Lagrangian)  bins  described  in  Section  2.4  is  employed.  In  this  case, 
the  bins  are  fixed  in  the  (initial)  reference  configuration  of  the  com¬ 
putational  cell.  In  other  words,  the  same  atoms  are  associated  with 
a  given  bin  throughout  the  entire  molecular  dynamics  simulation. 


Fig.  5.  (a)  Temporal  evolution  of  the  particle  velocity  associated  with  the  propaga¬ 
tion  of  two  approaching  shock  waves.  The  results  in  (b)  are  associated  with  a  50% 
higher  computational-cell  axial  contraction  rate. 

Examples  of  the  typical  results  obtained  through  the  use  of  this 
method  are  displayed  in  Fig.  5(a)  and  (b).  The  results  displayed  in 
Fig.  5(a)  and  (b)  are  obtained  under  identical  conditions  except  for 
the  rate  of  axial  contraction  of  the  computational  cell  (50%  higher 
in  the  case  of  Fig.  5(b)).  In  these  figures,  particle  velocities  at  differ¬ 
ent  simulation  (i.e.  post  shock  wave  generation)  times  are  plotted 
against  the  Lagrangian  bin  center  x-location.  Brief  examination  of 
the  results  displayed  in  these  figures  reveals  that: 

(a)  two  shock  waves  are  generated  (only  the  right-propagating 
shock  is  shown  though)  at  the  computational  cell  faces  nor¬ 
mal  to  the  x-direction.  These  shocks  then  propagate  towards 
the  computational-cell  center; 

(b)  after  a  brief  transient  period,  the  shocks  appear  to  reach  a  steady 
wave  profile  (i.e.  a  time-invariant  profile  within  a  reference 
frame  which  is  attached  to,  and  moves  with,  the  shock  front); 
and 

(c)  both  the  particle  velocity  and  the  shock  speed  increase  with 
the  computational-cell  contraction  rate.  It  should  be  noted  that 
the  curves  bearing  the  same  numerical  label  in  Fig.  5(a)  and  (b) 
correspond  to  the  same  simulation  time. 
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Fig.  6.  The  Hugoniot  relation  pertaining  to  the  particle  velocity  dependence  of  the 
shock  speed. 


It  should  also  be  noted  that  no  thermostat  was  used  in 
the  present  non-equilibrium  molecular  dynamics  simulations,  so 
that  the  steady-wave  shock  profile  is  a  natural  consequence 
of  a  balance  between  the  continuous  supply  of  momentum  to 
the  system  (through  the  continuous  computational  cell  axial 
contraction)  and  the  observed  lateral  motion  of  the  atoms  in 
the  continuously  enlarged  upstream  material  domain  swept  by 
the  shock.  In  general,  the  use  of  a  thermostat  modifies  the 
(F=ma  Newtonian-type)  equations  of  motion  solved  within  the 
molecular  dynamics  simulations  by  the  introduction  of  a  (velocity- 
proportional)  viscous-dissipation  term.  It  is  well-established  that, 
when  shock  formation  and  propagation  is  analyzed  within  a  contin¬ 
uum  framework,  the  use  of  a  viscous-dissipation  term  is  mandatory 
for  the  attainment  of  a  steady-wave  shock  profile.  This  fact  has 
often  been  used  as  a  justification  for  the  use  of  a  (local  or  global) 
thermostat  within  molecular-level  simulations  of  shock-wave  for¬ 
mation/propagation.  While  such  practices  greatly  facilitate  the 
attainment  of  a  steady  shock,  they  cannot  be  readily  defended  since 
shock  formation  and  propagation  is  generally  considered  to  be  an 
adiabatic  (no  system/surrounding  energy  exchange)  process  due  to 
the  near-instantaneous  material  transition  to  the  shocked  state.  In 
addition,  shock  generation/propagation  is  a  non-isentropic  process 
due  to  the  attendance  of  various  energy  dissipation  mechanisms. 
This  is  the  main  reason  that  no  thermostat  was  used  in  the  present 
work.  It  is  interesting  to  point  out  that  despite  the  fact  that  no 
viscous-dissipation  term  was  added  to  the  Newton’s  equations  of 
motion,  the  results  displayed  in  Fig.  5(a)  and  (b)  show  some  of  the 
defining  features  of  shock-waves  when  they  are  analyzed  in  the 
continuum-level  simulations  (in  the  presence  of  viscous  dissipa¬ 
tion).  Specifically,  a  steady  shock  is  generated  and  the  shock  width 
decreases  with  an  increase  in  the  shock  strength. 

3.2.  Determination  of  shock-Hugoniot  relations 

The  results  presented  in  Fig.  5(a)  and  (b)  reveal  the  steady-shock 
profile  and  can  be  used  to  find  a  functional  relation  between  the 
shock  speed,  Us ,  and  the  particle  velocity,  up.  This  was  done  in  the 
present  work  and  the  results  obtained  are  displayed  in  Fig.  6.  The 
two  curves  displayed  in  Fig.  6  labeled  respectively  as  “mixed”  and 
“segregated”  correspond  to  the  fully  mixed  and  micro-phase  segre¬ 
gated  initial  states  of  the  polyurea.  The  Us  vs.  up  relation  is  one  of  the 
so-called  shock  Hugoniot  relations.  In  general,  the  Hugoniot  can  be 
defined  as  a  locus  of  shocked-material  states  in  a  stress/pressure, 


energy  density,  mass  density  (or  specific  volume),  temperature, 
particle  velocity  and  shock  speed  space  which  are  associated  with 
(or  “centered  at”)  a  given  initial/reference  material  state.  The  Us  vs. 
up  relation  mentioned  above  is  a  simple  projection  of  the  Hugo¬ 
niot  to  the  Us-up  plane.  In  the  case  of  planar  shocks,  of  interest 
in  the  present  work,  the  other  commonly  used  Hugoniot  relations 
include:  axial  stress,  tn,  vs.  density,  p  (or  specific  volume,  v  =  1/p); 
(mass-based)  internal  energy  density,  e,  vs.  p  (or  v);  tn  vs.  up  and 
temperature  T  vs.  p  (or  v).  These  relations  were  determined  in  the 
present  work  using  two  distinct  methods: 


(a)  The  first  method  is  based  on  the  three  so-called  “jump  equa¬ 
tions”  which  are  defined  as: 


P~US  =  P+(US  -  Up)  (6) 

t-  +  p~Uf  =  t+  +  p+[Us  -  up)2  (7) 

t  £+ 

e~  +  —  +  0.5Uj  =e+  +  -TT  +  0.5(US  -  u„)2  (8) 

p~  p+ 

These  equations  relate  the  known  downstream  material 
states  (denoted  by  a  superscript  “-”)  and  the  unknown 
upstream  material  states  (denoted  by  a  superscript  “+”)  asso¬ 
ciated  with  the  shock  of  a  given  strength  (as  quantified  by  the 
shock  speed  or  the  downstream-to-upstream  particle  velocity 
jump).  These  equations  are  next  combined  with  the  previously 
determined  Us  vs.  up  relation  and  the  prescribed  (shock- 
strength  defining  quantity)  Us  or  up  to  solve  for  the  unknown 
upstream  material  states.  It  should  be  noted  that  this  method 
enables  determination  of  only  material  mechanical  state  vari¬ 
ables  (tn,  e,  v(=  1/p),  Us  and  up).  To  obtain  temperature,  a 
separate  set  of  equilibrium  NVE  (E-total  energy  of  the  sys¬ 
tem)  molecular  dynamics  simulations  is  carried  out.  In  each 
case,  a  local  computational  sub-cell  is  defined  containing  only 
the  upstream  (shocked)  material.  The  number  of  particles,  the 
volume  of  the  sub-cell  and  its  total  internal  energy  are  all  main¬ 
tained  constant.  The  associated  “ equilibrium ”  temperature  is 
then  calculated  using  the  time-averages  of  the  atomic  velocities 
(see  Eq.  (10)  below);  and 

(b)  Time  averages  of  the  atomic  positions,  rz,  velocities,  vit  and 
interaction  forces,/;  (i  is  the  atomic  label)  are  used  to  compute 
the  unknown,  local  (bin-based)  thermo-mechanical  quantities 
using  the  following  standard  thermodynamic  relations: 


P  = 


T  = 


Vbir 


avg 


E 

1  1  Ef?w'-"1 


L  i=i 


3Nbinkb 

Hi  =  77—  |  NbinkbT  +  V rf  ®fi 

Vbin  v  tr  , 

Z=Uotai  +  \(^rniV,vX) 


avg 


(9) 


(10) 


(11) 


(12) 


k  i=l 


avg 


where  subscript  “avg”  denotes  time  averaging,  Nbin  and  Vbin  the 
number  of  atoms  within  and  the  volume  of  the  bin,  respectively, 
kb  is  the  Boltzmann’s  constant,  ETotai  is  given  by  Eq.  (1 ),  while 
and  “(g)”  indicate  dot  product  and  tensorial  product  operators, 
respectively. 

It  should  be  noted  that,  in  this  case,  the  bins  were  defined 
within  a  reference  frame  which  is  attached  to  (and  moves  with) 
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the  steady-shock  front.  Clearly,  in  this  case  different  atoms  may 
reside  within  a  given  bin  at  different  simulation  times.  On  the 
other  hand,  the  bins  correctly  collect  the  information  about  the 
atoms  (temporarily)  residing  in  a  given  portion  of  the  steady- 
shock  profile.  In  other  words,  time  averages  are  calculated  not 
for  a  fixed  assembly  of  atoms,  but  rather  for  a  transient  set  of 
atoms  associated  with  a  given  moving  bin.  It  should  be  also 
noted  that  since  one  of  the  main  objectives  of  the  present  work 
was  determination  of  the  Hugoniot  relations,  only  the  data 
pertaining  to  the  bins  located  in  the  upstream  portion  of  the 
shock- wave  are  collected  and  analyzed  (for  different  shock- 
strength  conditions). 

The  results  of  these  two  procedures  are  displayed  in 
Fig.  7(a)-(d).  In  each  of  these  figures,  two  cases  are  shown  and 
labeled  as  “Method  ( a )”  or  “ Method  (b)”  in  order  to  indicate  the 
method  used  for  generation  of  the  corresponding  results.  In 
addition,  separate  curves  are  given  for  the  two  initial  conditions 
of  polyurea.  It  is  apparent  that  the  two  methods  yield  consistent 
results  while  the  initial  (“mixed”  or  “ segregated ”)  condition  of 
the  polyurea  has  a  noticeable  effect  on  to  the  computed  shock 
Hugoniot  relations.  An  explanation  for  the  latter  observation 
will  be  provided  in  the  next  section. 

The  Hugoniot  relations  displayed  in  Figs.  6  and  7(a)-(d)  are  typ¬ 
ically  used  within  a  continuum-level  computational  analyses  of 
shock- wave  generation/propagation  in  two  ways: 

(a)  They  are  directly  used  in  the  analyses  of  shock-wave  propaga¬ 
tion  under  uniaxial  strain  conditions  [e.g.  3];  and 

(b)  Alternatively,  they  can  be  used  to  derive  a  continuum-level 
material  model  which  is  consistent  with  the  material  mechan¬ 
ical  response  under  high-rate,  large-strain,  high-pressure 
conditions.  Such  a  model  is  subsequently  used  in  general  three- 
dimensional,  non-linear  dynamics  computational  analyses  [4]. 

In  our  recent  work  [20],  it  was  shown  that  the  Hugoniot  rela¬ 
tions  can  be  generated  by  converting  the  corresponding  isotherms 
(obtained  via  quasi-static,  molecular-level  mechanical  tests).  This 
procedure  was  found  to  be  associated  with  a  number  of  challenges 
(e.g.  a  particular  form  of  the  equation  of  state  had  to  be  assumed, 
several  material  properties/relations  had  to  be  assessed  indepen¬ 
dently,  etc.).  Most  of  these  challenges  were  not  encountered  in  the 
present  work  since  the  Hugoniot  relations  are  derived  more  directly 
from  the  molecular-level  computational  results. 

3.3.  Shock-induced  material-state  changes 

The  results  presented  and  discussed  in  the  previous  sections 
clearly  revealed  the  formation  and  propagation  of  planar  shocks 
in  polyurea  and  enabled  formulation  of  the  appropriate  shock- 
Hugoniot  relations.  In  addition,  a  clear  difference  in  the  shock- wave 
form  and  the  associated  Hugoniot  relations  was  observed  in  the 
cases  of  fully-mixed  and  micro-phase  segregated  polyurea.  In  the 
present  section,  a  more  detailed  investigation  is  carried  out  of  the 
molecular-level  material  microstructure  in  the  wake  of  a  propa¬ 
gating  planar  shock  in  order  to  provide  some  justification  for  the 
observed  differences  in  a  shock- wave  behavior  in  the  two  polyurea 
initial  states. 

It  should  be  first  realized  that  the  analysis  of  material 
microstructure  and  its  evolution  due  to  shock  loading  in  polyurea 
is  quite  challenging  due  to  the  absence  of  a  crystal  structure  in 
this  material.  Namely,  when  molecular-level  simulations  of  shock 
generation/propagation  are  carried  out  in  (nearly  perfect)  single 
crystal  solids  (e.g.  [8,11]),  deviations  from  long  range  order  (i.e. 
formation  of  various  point,  line  and  planar  defects)  can  be  readily 
investigated.  Polyurea  is,  on  the  other  hand,  an  amorphous  material 


in  its  initial  condition  and  remains  so  after  being  subjected  to  the 
shock  loading.  To  address  the  challenge  of  material  microstructure 
characterization  and  its  changes  resulting  from  shock  loading,  two 
microstructural  parameters  are  monitored  in  the  present  work:  (a) 
the  extent  of  lateral  motion  of  the  atoms  at  the  shock  front;  and  (b) 
the  extent  of  hydrogen-bond  breaking  (in  the  case  of  micro-phase 
segregated  polyurea,  only). 

In  the  case  of  single-crystalline  solids,  previous  molecular-level 
shock-wave  formation/calculation  work  (e.g.  [8,11])  established 
that  the  steady-wave  condition  is  attained  not  as  a  result  of  viscous 
dissipation  (as  is  the  case  for  shocks  in  fluids)  but  rather  as  a  result 
of  transverse  atomic  motions  which  result  in  inelastic  deformation 
(permanent  slippage  of  crystal  planes  and  the  formation  of  crystal 
defects).  Similar  lateral  motion  of  the  atoms  at,  or  in  the  vicinity 
of,  the  shock  front  is  also  observed  in  the  present  work.  Since  these 
motions  resulted  in  a  permanent  change  in  atomic  nearest  neighbor 
coordination,  they  can  be  considered  as  in-elastic  strain  producing. 
However,  as  stated  earlier,  the  amorphous  nature  of  polyurea  pre¬ 
cludes  a  more  detailed/quantitative  description  of  the  nature  of  the 
atomic  lateral  motion  process  or  the  associated  defect  formation. 
What  is  certain,  however,  is  that  no  covalent  bond-breaking  was 
observed  to  take  place.  In  addition,  shock-loading  is  found  to  lead 
to  a  permanent  (1-2%)  densification  of  this  material  (for  both  the 
fully-mixed  and  micro-phase  segregated  initial  conditions).  This 
finding  is  somewhat  puzzling  since  polyurea,  within  the  continuum 
framework,  is  considered  as  a  nearly  incompressible  material. 

Significant  lateral  atomic  motions  are  observed  in  both  the  fully 
mixed  and  micro-phase  segregated  states  of  polyurea.  No  statisti¬ 
cally  significant  differences  in  the  extents  of  this  motion  could  be 
detected  for  the  two  states  of  polyurea.  However,  numerous  obser¬ 
vations  of  hydrogen-bond  breaking  between  urea  linkages  of  the 
adjacent  chains  are  made  in  the  case  of  the  micro-phase  segre¬ 
gated  polyurea.  Since  hydrogen-bond  breaking  is  a  fairly  potent 
shock-wave  energy  absorbing/dissipating  mechanism,  one  may 
expect  that  this  process  contributes  significantly  to  the  experimen¬ 
tally  observed  high  blast/ballistic  impact  mitigation  potential  of 
polyurea. 

To  demonstrate  that  hydrogen-bond  breaking  leads  to  irre¬ 
versible  changes  in  the  material  properties,  molecular  statics-based 
local  simple-shear  tests  of  the  upstream  material  were  carried  out. 
Within  these  tests,  a  local  (upstream)  computational  cell  is  sub¬ 
jected  to  a  sequence  of  simple-shear  deformation  modes  and,  in 
each  case,  an  energy  minimization  procedure  is  employed.  Then, 
the  procedure  used  in  our  previous  work  [21-24],  which  relates 
the  (minimum)  energy  to  the  first  and  second  invariants  of  the 
left  Cauchy  deformation  tensor,  is  employed  in  order  to  determine 
the  evolution  of  the  stress  state.  The  initial  rate  of  change  of  the 
resulting  shear  stress  with  the  change  in  the  shear  strain  is  used 
as  a  measure  of  the  material  shear  modulus.  This  procedure  clearly 
revealed  that  hydrogen-bond  breaking  due  to  material  exposure 
to  shock  loading  leads  to  a  10-15%  reduction  in  the  material  shear 
modulus.  An  example  of  the  relaxed  molecular-level  microstruc¬ 
ture  of  micro-phase  segregated  and  shocked  polyurea  after  being 
subjected  to  simple  shear  is  displayed  in  Fig.  8. 

3.4.  Molecular-level  impulse  test 

To  help  reveal  differences  in  the  mechanical  response  of  fully 
mixed  and  micro-phase  segregated  polyurea  to  shock  loading,  a 
molecular-level  impulse  (computational)  test  is  carried  out.  Within 
this  test,  the  computational-cell  axial  lattice  parameter  is  perturbed 
only  once  by  a  finite  amount.  This  created  a  pair  of  impulse  waves 
(one  at  each  y-z  face  of  the  cell).  The  impulse  essentially  consists  of 
a  shock- wave  closely  followed  by  a  release  wave.  An  example  of  the 
typical  results  obtained  in  this  portion  of  the  work  is  given  in  Fig.  9 
in  which  impulse  profiles  at  different  simulation  times  are  shown 
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Fig.  7.  (a)  Axial  stress  vs.  specific  volume;  (b)  energy  vs.  specific  volume;  (c)  axial  stress  vs.  particle  velocity;  and  (d)  temperature  vs.  specific  volume  Hugoniot  relations. 
Please  see  text  for  explanation  of  the  “Method  (a)”  and  “Method  (b)”. 


Bin  Center  Location 


Fig.  8.  Relaxed  atomic  configurations  of  the  micro-phase  segregated  and  shocked 
polyurea  after  being  subjected  to  simple  shear. 


Fig.  9.  An  example  of  the  typical  results  obtained  in  the  present  work  pertaining  to 
the  temporal  evolution  (dispersion)  of  an  impulse  travelling  to  the  right. 
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for  the  right-travelling  impulse,  only.  As  is  seen  in  Fig.  9,  the  impulse 
tends  to  become  dispersed  as  it  advances  towards  the  center  of  the 
computational  cell.  Ultimately,  the  impulse  becomes  completely 
dispersed  (i.e.  difficult  to  resolve).  The  main  task  within  the  present 
section  was  to  determine  if  there  is  a  statistical  difference  in  the 
distance  of  propagation  of  the  impulses  in  the  two  polyurea  cases 
before  the  impulses  become  completely  dispersed. 

The  results  obtained  in  the  present  work  revealed  that  the 
micro-phase  segregated  polyurea  consistently  yields  a  shorter 
impulse  “ complete-dispersion ”  distance  (by  5-10%  depending  on  the 
shock  strength).  While  this  finding  is  quite  preliminary  and  requires 
further  in  depth  investigation,  it  provides  a  tentative  proof  that 
shock-induced  hydrogen-bond  breaking  plays  an  important  role  in 
the  superior  shock-impact  mitigation  capacity  of  this  material. 

3.5.  Final  remarks 

Within  the  present  work,  molecular-level  computational  meth¬ 
ods  are  employed  to  study  various  phenomena  accompanying 
shock-wave  generation  and  propagation  in  polyurea.  It  should  be 
noted  that  the  present  work  does  not  suggest  that  molecular- 
level  analyses  of  shock  generation  and  propagation  should  replace 
the  corresponding  continuum  (hydrodynamic)  analyses.  The  lat¬ 
ter  are  far  better  suited  (and  feasible)  for  studying  the  behavior 
of  real-life  engineering  systems  (e.g.  suspension  layers  within  the 
blast-protection  helmets).  Rather,  the  present  approach  is  highly 
beneficial  relative  to  the  identification  and  characterization  of  the 
nano-scale  phenomena/processes  taking  place  at  the  shock  front. 
Once  these  phenomena/processes  are  well  understood  and  char¬ 
acterized  (a  formidable  task),  the  knowledge  gained  can  be  used 
to  formulate  (and  parameterize)  more  physically  based  material 
models  suitable  for  use  in  continuum-level  analyses. 

Another  important  aspect  of  the  present  work,  worth  discussing, 
is  the  size  of  the  computational  cell  used.  This  size  corresponds  to 
the  maximum  size  that  can  be  handled  (within  realistic  computa¬ 
tional  times)  with  the  given  combination  of  the  number  of  software 
(Discover)  licenses  and  the  computer  resources  available  during  the 
present  work.  The  main  consequence  of  the  relatively  small  com¬ 
putational  cell  used  is  that  the  present  study  was  limited  to  the 
intermediate-  and  high-strength  shock  regimes  (i.e.  the  regimes 
which  are  associated  with  a  reasonably  small  shock  width).  These 
regimes  may  not  be  of  main  interest  within  the  context  of  traumatic 
brain  injury  since  they  would  (most  likely)  cause  soldier  fatality. 
Weaker  shocks,  which  are  not  too  strong  to  cause  a  soldier’s  death 
yet  strong  enough  to  cause  TBI,  are  of  major  concern.  These  shocks 
are  associated  with  a  much  larger  width  which  could  not  be  accom¬ 
modated  with  the  present  size  of  the  computational  cell.  We  are 
in  the  process  of  acquiring  computational  capabilities  which  will 
allow  us  to  address  the  weak  shock-wave  regime  in  the  future.  This 
will  allow  us  to  generate  the  appropriate  Hugoniots  and  identify 
the  nature  of  the  energy  dissipation  processes  in  this  shock-wave 
regime. 

4.  Summary  and  conclusions 

Based  on  the  results  obtained  in  the  present  work,  the  following 
summary  remarks  and  main  conclusions  can  be  drawn: 

Various  phenomena  accompanying  the  formation  and  prop¬ 
agation  of  a  planar  shock-wave  within  polyurea,  a  micro-phase 
segregated  elastomer,  are  investigated  using  molecular-level  com¬ 
putational  methods. 

The  results  obtained  show  that  even  without  the  use  of  a 
viscous-dissipation-based  thermostat,  a  steady-wave  planar  shock 
profile  can  readily  be  established  in  this  material. 

The  time-averaged  results  pertaining  to  the  atomic  positions, 
velocities  and  interaction  forces  are  used  to  construct  the  appropri¬ 


ate  shock  Hugoniot  relations,  the  relations  which  define  the  locus 
of  stress,  energy,  density,  temperature  and  particle  velocity  of  the 
material  swept  by  a  shock  propagating  at  a  given  speed. 

Detailed  examination  of  the  molecular-level  microstructure 
evolution  in  the  shock-wave  wake  is  carried  out  in  order  to  identify 
the  nature  of  energy  absorbing  and  shock-wave  spreading  mech¬ 
anisms.  The  results  revealed  that  shock  loading  causes  extensive 
hydrogen  bond  breaking  in  the  micro-phase  segregated  polyurea. 
These  processes  are  associated  with  substantial  energy  absorption 
and  dissipation  and  are  believed  to  be  related  to  the  experimen¬ 
tally  observed  high  blast/ballistic  impact  mitigation  potential  of 
polyurea. 
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